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ABSTRACT 

This  paper  proposes  a  hierarchical  Bayesian  model  for  multiple-pass,  multiple  antenna  synthetic  aperture 
radar  (SAR)  systems  with  the  goal  of  adaptive  change  detection.  We  model  the  SAR  phenomenology  directly, 
including  antenna  and  spatial  dependencies,  speckle  and  specular  noise,  and  stationary  clutter.  We  extend 
previous  work1  by  estimating  the  antenna  covariance  matrix  directly,  leading  to  improved  performance  in 
high  clutter  regions.  The  proposed  SAR  model  is  also  shown  to  be  easily  generalizable  when  additional  prior 
information  is  available,  such  as  locations  of  roads/intersections  or  smoothness  priors  on  the  target  motion. 
The  performance  of  our  posterior  inference  algorithm  is  analyzed  over  a  large  set  of  measured  SAR  imagery. 
It  is  shown  that  the  proposed  algorithm  provides  competitive  or  better  results  to  common  change  detection 
algorithms  with  additional  benefits  such  as  few  tuning  parameters  and  a  characterization  of  the  posterior 
distribution. 
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1.  INTRODUCTION 

SAR  radars  traditionally  are  used  to  image  stationary  objects  by  integrating  radar  returns  from  spatially 
diverse  points.  SAR  images  are  often  formed  with  much  longer  integration  times  than  other  radar  modes  (in 
particular,  MTI-mode  radars).  While  this  leads  to  high  resolution  images  in  many  cases,  the  long  integration 
times  also  lead  to  other  undesirable  effects.  Moving  objects  can  cause  phase  errors  in  the  reconstruction  of 
SAR  images,  leading  to  well  known  dispersion  and  displacement  effects.2,3  Moreover,  angular  scintillation 
(aka  ‘specular  noise’  or  ’glints’)  can  cause  large  phase  errors  that  may  considerably  degrade  signal  quality.4 
As  the  name  suggests,  angular  scintillation  has  a  large  angular  dependence  in  the  sense  that  the  intensity  of 
this  noise  source  is  only  large  from  few  azimuth  angles.  Speckle  noise  is  an  additional  noise  source  that  arises 
from  coherent  imaging  with  SAR.  Speckle  noise  tends  to  be  spatially  correlated,  depending  on  the  texture 
of  the  surrounding  pixels  (e.g.,  buildings  versus  vegetation.)  Unlike  specular  noise,  the  intensity  of  speckle 
noise  tends  to  be  uniform  as  a  function  of  azimuth  angle. 

In  this  work,  we  focus  on  the  situation  where  the  scene  is  being  persistently  monitored.  Thus,  images  are 
available  from  multiple  passes  of  the  radar  platform,  multiple  antennas  (phase  centers),  and  multiple  look- 
angles.  However,  the  background  of  these  images  remains  relatively  unchanged  and  is  thus  modeled  as  being 
embedded  in  a  low-dimensional  subspace.  Exploiting  the  low-dimensional  subspace  in  order  to  extract  the 
moving  targets  of  interest  is  used  by  many  common  algorithms.  Soumekh5  shows  that  all  stationary  objects 
can  be  removed  from  an  ideal  set  of  SAR  images  formed  in  a  monopulse  radar  system  by  using  a  simple 
difference  image,  also  known  as  displaced  phase  center  array  (DPCA)  processing.  In  practice,  the  antennas 
are  not  perfectly  calibrated,  which  can  significantly  degrade  the  performance  of  DPCA.  Signal  subspace 
processing6  (SSP)  addresses  this  problem  by  proposing  an  adaptive  blind  calibration  technique  of  a  two 
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channel  system  that  assumes  that  the  signals  are  related  through  a  locally  spatially  invariant  filter.  Ender 
also  considers  using  multiple  channels  for  SAR  detection  of  moving  targets  through  space-time  adaptive 
processing'  (STAP),  which  adaptively  estimates  the  clutter  covariance  matrix  in  order  to  separate  moving 
targets  from  the  background.  However,  STAP  is  limited  in  the  sense  that  it  requires  target-free  data  in  order 
to  effectively  estimate  the  clutter  covariance  matrix.  In  our  work,  we  estimate  both  the  calibration  constants 
and  the  clutter  covariance  directly  in  order  to  extract  the  moving  targets  from  the  SAR  images. 

Previously,  we  presented  a  model  for  synthetic  aperture  radar  (SAR)  imagery  in  order  to  separate  a  sparse 
foreground  component  that  may  contain  targets  of  interest  from  the  stationary  background.1  A  hierarchical 
Bayesian  formulation  was  proposed  that  extends  the  model  of  Ding,  He  and  Carin8  using  the  methods  of 
Tipping9  and  Wipf  and  Rao. 10,11  This  model 

1.  Provided  Monte  Carlo  estimates  of  the  posterior  distribution  using  an  efficient  MCMC  sampler. 

2.  Used  conjugate  distributions  with  non-informative  hyperparameters,  leading  to  very  few  tuning  param¬ 
eters. 

3.  Had  low  reconstruction  errors  and  robustness  to  noise  parameters  on  a  synthetic  dataset. 

4.  Performed  well  on  a  limited  set  of  measured  SAR  imagery. 

In  this  paper,  we  extend  this  model  in  several  ways  that  include 

1.  Estimating  the  antenna  covariance  matrix  directly,  leading  to  improved  performance  in  both  low  and 
high  variance  regions  of  the  scene. 

2.  Generalizing  the  model  to  be  able  to  account  for  prior  knowledge  on  the  locations  of  the  target  signature. 

3.  Generalizing  the  model  to  include  spatial  dependencies  using  a  hidden  Markov  model. 

4.  Analyzing  the  performance  of  our  algorithm  over  a  large  set  of  measured  SAR  imagery  with 

•  Comparisons  to  common  change  detection  algorithms  such  as  DPCA  and  along-track  interferom¬ 
etry  (ATI). 

•  Comparisons  as  a  function  of  our  model  (i.e. ,  with/without  prior  knowledge  and/or  an  antenna 
covariance  model). 

•  Analysis  as  a  function  of  integration  length  and  the  number  of  observations. 

The  rest  of  the  paper  is  organized  as  follows.  The  SAR  image  model  is  given  in  Section  2.  Section  3 
discusses  the  Monte  Carlo  implementation  for  posterior  inference  on  the  model.  Performance  of  the  model 
and  inference  is  analyzed  with  both  synthetic  and  measured  SAR  imagery  in  Section  4.  Finally,  conclusions 
and  future  work  are  provided  in  Section  5. 

2.  SAR  IMAGE  MODEL 

In  this  work,  it  is  assumed  that  we  have  access  to  SAR  images  formed  from  multiple  passes  of  a  radar  platform 
with  multiple  antennas  (i.e.,  phase  centers.)  Moreover,  images  are  formed  over  distinct  azimuth  angle  ranges 
that  can  be  indexed  by  the  frame  number,  /.  Table  1  provides  the  indexing  scheme  used  throughout  this 
paper  in  order  to  distinguish  between  images  from  various  antennas,  frames,  and/or  passes. 

The  image  model  proposed  in  this  paper  is  very  similar  to  our  previous  work.1  Whenever  relevant,  we 
will  point  out  the  differences  in  our  models. 


Table  1.  Index  variable  descriptions 


Index  Description 

Index  Variable 

Range 

Antenna  (channel) 

k 

1,2,.. 

.,K 

Frame  (azimuth  range) 

f 

1,2,. 

■,F 

Pass 

i 

1,2,.. 

;N 

Pixel 

V 

1,2,. 

,.,P 

(c)  Locally  identical  distributions  (large)  (d)  HMM  dependence  for  a  single  pixel 

Figure  1.  Various  scenarios  for  spatial  dependence  used  in  this  paper.  In  (a),  pixels  are  independent.  In  (b)  and  (c), 
pixels  share  an  identical  distribution  over  a  medium  and  large  area,  respectively.  In  (d),  pixels  are  assumed  to  have 
a  hidden  Markov  model  dependence  on  neighboring  cells  with  higher  probabilities  close  to  the  center  pixel. 


Similar  to  our  previous  work,  we  propose  a  decomposition  of  SAR  images  into  a  background  component 
that  lies  in  a  low-dimensional  subspace,  a  foreground  component  that  is  sparse,  and  additive  noise: 

=  HkjAP)  ( Bk,fAP )  +  Tk,f,i(p ))  +  Vkj,i(p),  (1) 

where  Hkj,i(p)  is  a  spatially  varying  filter,  Bkj,i{p)  is  the  background  component,  Tk.jhp)  is  the  foreground 
component,  and  14 J,i(p)  is  zero-mean  additive  noise.  Moreover,  we  decompose  the  background  component 
as 

Bk,f,i(p)  =  skj(p )  +  Xk,f,i(p)  +  Sf(p)Yk,fli(p),  (2) 

where  Skj{p)  is  the  inherent  background  that  is  identical  over  all  passes,  Xkjj{p)  is  the  speckle  noise 
component,  Y^ja ip)  is  the  specular  noise  (glint)  component,  and  5f(p)  is  an  indicator  variable  for  specular 
noise.  The  foreground  component  decomposes  as 

Tk,f,i(p)  =  Dfti(p)Mkifti(p),  (3) 

where  is  the  target  value  and  Df  i(p)  is  an  indicator  variable  for  the  existence  of  the  target  at  that 

pixel.  Table  2  provides  the  distributional  models  of  each  of  these  variables.  The  remainder  of  the  section 
discusses  the  model  in  greater  detail. 


Table  2.  Distributional  models  for  each  component  in  equations  (1),  (2),  and  (3) 


Component 

Variable 

Distribution  | 

Parameters 

Region 

All  jfc? 

Calibration  constant 

Deterministic  (LS) 

none 

PGZg 

No 

Stationary 

Skj{p) 

Normal  (0,  Es(j)) 

£s(j) 

p  e  Qj 

Yes 

Speckle 

Xk,f,i(p) 

Normal  (0,  £A  (j)) 

xxU) 

p  G  Qj 

Yes 

Specular  (glint) 

Normal  (0,£^  (j)) 

sy(j) 

p  G  Qj 

Yes 

Specular  indicator 

hip) 

Bernoulli  (^ttJ  (p)'j 

7T/  (?) 

Each  p 

No 

Target 

Mk,f,i(p) 

Normal  (0,  £M) 

SM 

All  p 

Yes 

Target  indicator 

Df,i{P ) 

Bernoulli  ^7 rj^(p)^ 

*%(P) 

Each  p 

No 

Additive  noise 

VkJ,i(p) 

Normal  (0,  Sv) 

Ev 

All  p 

Yes 

Class  assignment 

dip) 

Multinomial(l;  q) 

Q 

Each  p 

No 

2.1  Calibration  constants,  Hkjj{p) 

The  calibration  constants  are  assumed  to  be  constant  within  small  spatial  regions  as  denoted  by  either 
Figure  1(b)  or  1(c).  The  calibration  constant  within  such  a  region  is  a  single  value  /ifc,/y(g)  for  all  pixels 
p  G  Zg.  In  general,  hkj,i(g)  may  be  randomly  distributed,  though  experimental  evidence  suggests  that  using 
least-squares  is  sufficient  for  estimating  the  values. 

2.2  Antenna  covariance  models  for  £ 

Each  pixel  in  the  stationary,  speckle,  specular,  and  target  value  components  is  jointly  distributed  across 
the  K  antennas.  Since  the  pixel  values  are  complex,  the  joint  distribution  is  over  2 K  values.  Consider  a 
zero-mean  vector  W  =  [u>i  W2  •  •  ■  wk]T,  where 


Wk  =  ak+jbk,  k=l,2,...,K. 


(4) 


In  this  work,  W  could  be  representative  of  {Skj(p)}k,  {Xkj,i(p)}k,  {Ykj,i{p)}k  or  {Mkjti(p)}k.  The  complex 
covariance  matrix  is  given  by 

r  =  E  [WWH]  €  CKxK,  (5) 

which  can  be  related  to  the  real- valued  covariance  matrix 


Re{T}  —Im  {T } 

/to  {r}  i?e{r} 


^  x2 K 


(6) 


In  this  work,  we  assume  a  specific  structure  for  the  complex  covariance  matrices  T  which  can  be  derived  by 
assuming  that 


E[(wk){wj)*] 


cr2pe  ^ kj 


k  =  j 


k,j  =  1,2, ...  ,K 


(7) 


and 


k,j  =  1,2, ...  ,K 


(8) 


where  is  a  rotation  matrix  by  angle  <f>.  In  particular,  let 


i 

pe-M 12  • 

■  perj<t>1K 

2 

pe-0<t>21 

1 

■  pe~^2K 

a 

pe~^K1 

pe~^K2  ■ 

■  pe-j4,KK 

where  a1  is  the  channel  variance,  p  is  the  coherence  between  antennas,  and  {< t>nm}n  m  are  the  phase  differ¬ 
ences  between  the  channels.  Note  that  for  background  components,  p  should  be  near  one  and  <f>ij  «  0.  A 
more  general  model  could  account  for  different  channel  variance  and  coherence  values,  but  since  we  use  the 
calibration  constants  Hk  f  j(p)  to  equalize  the  channels,  the  effect  was  seen  to  be  relatively  insignificant. 

It  should  be  noted  that  this  form  of  the  covariance  matrix  is  related  directly  to  some  of  the  very  common 
methods  for  change  detection  in  SAR  imagery.  In  particular,  consider  the  two  antenna  case  ( K  =  2).  It  can 
easily  be  shown  that  the  eigendecomposition  of  £  leads  to  eigenvalues  A  and  eigenvectors  v: 


A(S)  =  {2a2(l+p),2(r2(l-/9)} 


K£)  = 


(10) 

(11) 


When  the  phase  (f>  is  zero,  the  second  eigenvector  reduces  to  [1  —  1]T,  which  can  be  interpreted  as  DPCA 
for  a  two-antenna  system.  Moreover,  along-track  interferometry  (ATI),  which  thresholds  the  phase  <f>,  clearly 
depends  on  the  eigenvalues  in  a  direct  manner.  Deming12  shows  that  ATI  performs  well  when  canceling 
bright  clutter  (i.e.,  high  er2  and  p  «  1),  while  DPCA  performs  well  for  canceling  dim  clutter  (i.e.,  small  a2 
and  p  !=s  0.)  In  our  work,  we  hope  to  gain  the  discriminating  power  of  both  DPCA  and  ATI  by  modeling  the 
covariance  matrices  directly. 

Furthermore  it  should  be  noted  that  this  covariance  model  generalizes  our  previous  model1  which  set 
p  =  1  and  4>nm  =  0  for  all  n  ^  m  for  the  stationary,  speckle  and  specular  components.  In  Section  4,  we 
compare  the  performance  of  both  models  with  measured  SAR  imagery. 


2.3  Spatial  dependencies  for  the  background  components 

In  our  previous  model,1  we  assumed  that  objects  in  the  background  can  be  defined  by  one  of  C  classes  (e.g., 
roads,  vegetation,  or  buildings.)  Let 

d{p )  =  {dc}c=  i  ~  Multinomial(l;  qi,  cfe,  •  •  ■  >  Qc)  (12) 

where  qn  is  the  prior  probability  of  the  c-th  region  type.  Then  the  class  assignment  is  the  single  location  in 
d  with  value  equal  to  one.  Since  object  classes  tend  to  have  high  spatial  dependencies,  the  previous  model 
assumed  that  pixels  within  a  small  neighborhood  all  had  a  single  class  type.  In  this  work,  we  introduce  a 
hidden  Markov  model  dependency  (see  Figure  1(d))  which  assumes  that  pixels  in  neighboring  locations  have 
a  high  probability  of  being  the  same  class.  This  small  distinction  provides  two  advantages:  (a)  we  relax 
the  problem  of  tuning  the  size  of  the  spatial  dependency  region;  and  (b)  we  allow  finer  transitions  between 
classes  on  a  pixel-by-pixel  basis.  Future  work  plans  to  generalize  the  model  to  allow  for  mixtures  of  object 
classes. 

2.4  Priors  on  the  distribution  parameters 

In  general,  the  parameters  of  the  distributions  in  Table  2  are  not  known  a  priori.  Instead,  we  propose  a 
hierarchical  Bayesian  model  based  on  the  work  by  Tipping9  and  others.8, 10, 11  In  this  model,  we  assume 
that  the  parameters  of  the  distribution  are  also  random  variables  that  need  to  be  estimated  from  the  data. 
In  most  cases,  we  choose  a  non-informative  distributions  on  the  parameters  in  order  to  implement  inference 
algorithms  with  few  tuning  parameters.  Table  3  shows  the  distributions  of  the  so-called  ‘hyper-parameters’ 
used  in  this  work.  For  this  work,  we  assume  that  the  covariance  matrices  are  given  by  equations  (6)  and  (9) 
so  that  they  are  parameterized  by  a2,  p ,  and  {( t>nm}n  m-  Unlike  the  other  parameters,  using  this  structure 
suggests  an  informative  prior  based  on  our  knowledge  of  the  phenomenology  of  the  SAR  sensor.  Future  work 
will  compare  to  the  non-informative  case,  where  £  is  Inverse-Wishart  distributed. 

As  in  Tipping,9  we  assume  that  the  hyperparameters  on  the  Inverse  Gamma  distributions  are  non- 
informative  leading  to 


aa  =  bCT  =  1(T6.  (13) 

On  the  other  hand,  since  the  background  coherence  should  be  near  unity,  we  let 

ap  =  0.9, 
bp  =  0.1. 

It  is  also  assumed  that  all  object  classes  have  equal  prior  probability,  so  that 

c„  =  1/C,  n  =  1,  2, . . . ,  C  (16) 


(14) 

(15) 


2.4.1  Indicator  probability  models 


We  have  a  great  deal  of  flexibility  for  choosing  the  hyperparameters  for  the  indicator  probabilities.  To  ensure 
that  only  a  small  percentage  of  the  pixels  contain  either  specular  noise  or  targets,  one  should  set  the  Beta 
parameters  so  that 


a(p) 

a(p)  +  b(p) 


0. 


(17) 


However,  in  many  cases  there  is  a  high  degree  of  spatial  and  temporal  dependence  for  the  indicator  variables. 
For  example,  both  targets  and  glints  have  spatial  extents  that  spread  over  several  pixels,  while  moving  targets 
tend  to  transition  smoothly  to  neighboring  pixels  in  sequential  frames.  Furthermore,  it  may  be  possible  to 
identify  regions  of  the  scene  where  specular  noise  or  targets  have  a  high  likelihood,  e.g.  at  edges  of  buildings 


Table  3.  Distributional  models  for  parameters  of  distributions  in  Table  2 


iii  the  former  case  and  near  road  intersections  in  the  latter  case.  This  information  may  be  learned  adaptively 
through  observed  data  or  through  human-in-the  loop  processing.  In  any  of  these  cases,  one  would  expect 


that 

«0) 

a(p)  +  b(p) 


(18) 


Define  Wj  ( p,S )  and  W^(p,D)  to  be  functions  that  map  the  indicator  variables  6  and  D ,  respectively,  to 
a  real  number.  For  example,  this  may  be  the  average  number  of  non-zero  indicators  in  the  neighborhood  of 
pixel  p ,  or  it  may  be  a  weighted  version  similar  to  Figure  1(d).  Then  define: 


'ahigh,  Wj(p,Sf)  >  £]paUal 

bhigh  i 

Wj  ip,  Sf  )  >  Spatial 

«/  ip)  =  < 

,  b)(p)=< 

(19) 

^  & low  5  else 

^  blow ) 

else 

and 


M 


ahigh,  W%(p,  Df,)  >  Sfpatial  and  W%ip,  £>/-!,*)  >  temporal 


i(p)  = 


^ low  1 


(20) 


else 


KM  = 


bhzgh ,  (p,  Dfii)  >  £fpatial  and  Wft  ip,  >  £femporal 


(21) 


else. 

Note  that  ( aiow ,  biow )  should  be  chosen  to  satisfy  equation  (17)  and  (ahigh,  bhigh)  should  be  chosen  to  satisfy 
equation  (18). 


2.5  Generalizability  of  the  model 

The  authors  would  like  to  stress  that  the  model  presented  here  is  easily  generalized  to  include  additional 
information  sources.  For  example,  if  a  target’s  state  (e.g.  position  and  velocity),  then  one  could  predict  the 
location  of  its  signature  within  the  scene  with  high  probability  (see  Jao2  or  Newstadt  et.  al.13).  Moreover, 
one  might  consider  different  target  classes  or  explicitly  estimating  the  target  phase  which  is  proportional  to 
its  radial  velocity. 


3.  POSTERIOR  INFERENCE 

In  this  work,  the  computation  of  the  posterior  distribution  estimation  is  done  using  a  Gibb  sampling  scheme, 
where  the  vast  majority  of  variables  can  be  easily  sampled  from  simple  distributions  such  as  the  multivariate- 
normal  and  beta  distributions.  Since  we  wish  to  use  the  structure  of  the  antenna  covariance  given  in  equation 
(9),  sampling  is  more  difficult  and  requires  other  sampling  methods,  such  as  the  Metropolis-Hastings  algo¬ 
rithm.  Appendix  A  provides  pseudocode  for  the  implementation  of  our  sampling  method. 

4.  PERFORMANCE  ANALYSIS 

Figure  2  shows  the  full  scene  used  for  our  performance  analysis.  The  images  are  formed  from  a  superset  of 
the  data  available  through  the  public- released  SAR-GMTI  challenge  problem.14  Images  were  formed  for  a  3 
second  interval  from  three  passes  of  the  radar  platform  and  for  each  of  three  antennas.  We  considered  five 
separate  locations  that  contained  the  signature  of  the  Durango  truck  from  the  SAR-GMTI  challenge  problem 
at  times  tf,  f  =  1,2,. ..,5.  Figure  3  shows  images  of  the  scenes  created  at  each  of  the  times  tf.  These  images 
contain  a  variety  of  scenarios,  including  stationary  targets,  weak  target  signatures,  strong  target  signatures 
{(4,4)}  and  barely  visible  target  signatures. 


Figure  2.  The  full  scene  used  for  the  performance  analysis  is  shown  in  (a),  while  the  yellow  boxes  in  (b)  correspond 
to  the  individual  scenes  used  for  comparison  of  various  scenarios.  Note  that  the  images  formed  in  this  dataset  are  a 
superset  of  the  public-released  SAR-GMTI  challenge  problem  dataset.14 


Figure  4  shows  the  output  of  our  algorithm  for  the  scene/time  combinations,  {(1, 1),  (2,  2),  (3,  3),  (4, 4),  (5,  5)} 
(i.e. ,  the  scenes  containing  the  Durango  truck  signature),  as  well  as  comparisons  to  DPCA  and  a  mixture  of 
DPCA  and  ATI.  For  clarity,  a  grayscale  image  refers  to  a  deterministic  output  such  as  the  original  image  or 
DPCA  output,  while  a  colored  image  refers  to  a  probability  or  variance  parameter.  It  can  be  seen  that  our 
algorithm  performs  at  least  as  well  both  DPCA  and  the  DPCA/ATI  mixture,  although  usually  with  fewer 
false  alarms.  Moreover,  both  of  the  alternative  algorithms  required  tuning  of  the  threshold  parameters  (for 
each  of  the  5  scenes).  Finally,  our  method  provides  estimates  of  the  posterior  distribution,  such  as  the  vari¬ 
ance  of  the  background  (which  could  be  used  for  stronger  hypothesis  testing  procedures)  and  the  probability 
of  target  existence. 

Figure  5  compares  the  performance  of  our  algorithm  for  scenes  2  and  3  as  a  function  of  the  number 
of  observations.  In  particular,  images  were  formed  for  each  scene  from  n  distinct  look-angles,  where  n  G 
{10,20,30,40,50}.  Each  ‘observation’  consists  of  a  set  of  images  formed  for  K  =  3  antennas  and  N  =  3 
passes.  It  can  be  seen  that  as  the  number  of  observations  grows,  the  performance  of  our  algorithm  improves 
with  fewer  false  alarms. 

In  Figure  6,  we  compare  the  performance  for  scenes  2  and  3  as  a  function  of  the  integration  length 
({0.2, 0.5, 1.0}  seconds)  used  to  form  the  images.  It  should  be  noted  that  for  each  integration  length,  a  set 
of  images  was  formed  that  encompassed  the  entire  3  second  interval  (although  possibly  with  fewer/greater 
number  of  frames,  F).  Thus,  this  simulation  compares  the  robustness  of  our  algorithm  to  using  an  arbitrary 
integration  length.  The  results  suggest  that  the  performance  of  our  algorithm  is  invariant  to  the  choice  of 
integration  length.  This  relaxes  the  selection  of  another  tuning  parameter  that  can  affect  performance  in 
other  algorithms  such  as  DPCA  and  ATI. 
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Figure  3.  Image  matrix  of  the  scenes/times  evaluated  in  this  performance  analysis.  The  signature  of  a  single  target 
(Durango  Truck)  is  known  to  exist  along  the  diagonal,  although  other  targets  are  in  other  scenes.  Rows  of  the  matrix 
correspond  to  the  same  location  as  in  Figure  2(b),  while  columns  correspond  to  different  times  (i.e.,  azimuth  ranges). 
It  should  be  noted  that  there  are  stationary  targets  {(1,1)},  weak  target  signatures  {(2,  2),  (3, 3)},  strong  target 
signatures  {(4,4)}  and  barely  visible  target  signatures  {(5,5)}. 


Figure  4.  This  figure  shows  the  output  of  our  algorithm  for  the  scene/time  combinations  for  the  scenes  containing  the 
Durango  truck  signature,  {(1, 1),  (2,  2),  (3,  3),  (4, 4),  (5,  5)}.  Grayscale  images  refer  a  deterministic  output  such  as  the 
original  image  or  DPCA  output,  while  colored  images  refer  to  a  probability  or  variance  parameter.  It  is  seen  that 
our  algorithm  performs  at  least  as  well  both  DPCA  and  the  DPCA/ ATI  mixture,  although  usually  with  fewer  false 
alarms. 


Number  of  Observations 


Original  10  20  30 


40  50 


2 


0} 


OJ 

CJ 

m 


3 


Figure  5.  This  figure  compares  the  performance  of  our  proposed  method  for  scenes  2  and  3  as  a  function  of  the  number 
of  observations  ,  n.  A  single  observation  consists  of  a  set  of  images  formed  at  a  set  of  azimuth  (look-)  angles  9n  for  each 
of  K  =  3  antennas  and  N  =  3  passes.  The  number  of  observations  was  varied  for  n  G  {10,20,30,40,50}.  It  can  be 
seen  that  as  n  grows,  the  performance  of  our  algorithm  improves  with  fewer  false  alarms.  Moreover,  the  performance 
of  the  algorithm  is  relatively  consistent  for  larger  n,  which  suggests  that  the  estimation  of  the  background  distribution 
can  be  done  with  relatively  few  observations. 
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Figure  6.  This  figure  compares  the  performance  of  our  proposed  method  for  scenes  2  and  3  as  a  function  of  the 
integration  length  ({0.2,  0.5,  1.0}  seconds)  used  to  form  the  images.  The  results  suggest  that  the  performance  of  our 
algorithm  is  invariant  to  the  choice  of  integration  length,  which  relaxes  the  selection  of  another  tuning  parameter  that 
can  affect  performance  in  other  algorithms  such  as  DPCA  and  ATI. 


Figure  7  shows  the  output  of  our  algorithm  both  with  and  without  the  new  antenna  covariance  model.  It 
can  be  seen  that  using  the  covariance  model  significantly  reduces  the  number  of  false  alarms  in  areas  where 
the  background  has  high  variance. 

Figure  8  shows  the  output  of  our  algorithm  when  prior  information  on  the  location  of  the  targets  might 
be  available.  For  example,  in  scene  1,  targets  are  likely  to  be  stopped  at  an  intersection.  We  show  the 
performance  improvement  for  two  scenes  that  contain  targets  in  our  high  probability  region  (mission)  and 
that  don’t  contain  targets  in  those  regions  (reference).  It  can  be  seen  that  by  including  the  prior  information, 
we  are  able  to  detect  targets  with  much  better  accuracy,  including  the  stationary  targets  in  scene  1.  On  the 
other  hand,  we  don’t  have  significant  performance  decreases  in  the  reference  scenarios. 

5.  CONCLUSIONS  AND  FUTURE  WORK 

In  this  paper  we  extended  the  development  and  analysis  of  a  hierarchical  Bayesian  model  for  persistent 
SAR  imagery,  along  with  a  Gibbs  sampling  scheme  to  efficiently  estimate  the  posterior  distribution.  This 
algorithm  can  infer  statistics  of  the  noise  without  extensive  tuning  of  hyperparameters,  yet  also  provides  a 
characterization  of  its  uncertainty  through  a  posterior  distribution.  The  previous  model  was  extended  by 
estimating  the  antenna  covariance  matrix  directly,  leading  to  improved  performance  in  high  clutter  regions. 
The  model  is  also  shown  to  be  easily  generalizable  when  additional  prior  information  is  available,  such  as 
locations  of  roads/intersections  or  smoothness  priors  on  the  target  motion.  The  performance  of  our  posterior 
inference  algorithm  is  analyzed  over  a  large  set  of  measured  SAR  imagery.  It  is  shown  that  the  proposed 
algorithm  works  at  least  as  well  as  common  change  detection  algorithms.  Moreover,  the  algorithm  is  shown 
to  work  well  in  a  number  of  cases  that  include  detecting  stationary  targets  at  intersections,  detecting  targets 
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Figure  7.  This  figure  compares  the  performance  of  our  proposed  method  both  with  the  new  antenna  covariance  model 
and  as  well  the  simpler  model  where  the  coherence  of  the  channels  is  set  to  p  =  1.  It  can  be  seen  that  using  the  more 
general  covariance  model  significantly  reduces  the  number  of  false  alarms  in  areas  where  the  background  has  high 
variance  (i.e,  high  clutter  regions.) 
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Figure  8.  This  figure  compares  the  performance  of  our  proposed  method  with  and  without  priors  on  target  signature 
locations.  In  scene  1,  targets  are  likely  to  be  stopped  at  an  intersection,  while  in  scene  2,  the  target  signature  is 
assumed  to  lie  in  a  low  clutter  region.  We  show  the  performance  improvement  for  two  scenes  that  contain  targets  in 
the  high  probability  region  (mission  scenes)  and  for  two  scenes  that  don’t  contain  targets  in  those  regions  (reference 
scenes).  It  can  be  seen  that  by  including  the  prior  information,  we  are  able  to  detect  targets  with  much  better 
accuracy,  including  the  stationary  targets  in  scene  1.  On  the  other  hand,  we  don’t  have  significant  performance 
decreases  in  the  reference  scenarios. 


iii  both  bright  and  dim  clutter  situations,  and  working  well  regardless  of  the  integration  length  used  to  form 
the  SAR  images. 

Future  work  will  include  the  development  of  algorithms  that  exploit  the  use  of  a  posterior  distribution  for 
improved  performance  in  a  signal  processing  task,  e.g.  detection,  tracking  or  classification.  In  particular,  we 
are  interested  in  using  algorithms  for  simultaneously  detecting  and  estimating  targets  over  a  sparse  scene  with 
resource  constraints15, 16  as  well  determining  the  fundamental  performance  limits  of  a  SAR  target  tracking 
system.  Furthermore,  we  would  also  like  to  consider  other  generalizations  to  our  model,  such  as  complex 
target  maneuvers,  multiple  target  classes,  and  explicit  tracking  of  the  target  phase. 
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for  iteration  =  1  to  Nburn_ln  +  Nsampies  do 
I  =  I./H  %  Image  Calibration 

for  all  frames  /  do 
%  BACKGROUND 

(S,  X,  S,  Y)=  updateBackground(/,  S ,  X ,  <5,  Y.  M,  B ,  0,  /) 

%  FOREGROUND 

(D,  M)=  updateForeground(/,  S ,  X,  Y,  M,  R,  0,  /) 

%  CALIBRATION  CONSTANTS 

for  all  calibration  regions  g ,  antennas  fc,  passes  i  do 

A  =  [Sk,f  (p)  +  XkJA  (p)  +  Sf  (p) .*Yktfii  ( p)]PGZg ,  b  =  [Ikj,i  (p)]Pezg 
hk,f,i{g)  =  lscov(A,  b )  %  Least  squares  solution  to  Ah  =  b 
=  hk,f,i(g),  Vp&Zg 
end  for  %  (regions,  antennas,  passes) 
end  for  %  (frames) 

%  HYPERPARAMETERS 

0  =updateHyperParameters(/,  S,  X ,  <5,  Y,  M,  B,  H) 
end  for  %  (iterations) 

Figure  9.  Gibbs  Sampling  Pseudocode 
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APPENDIX  A.  PSEUDOCODE  FOR  POSTERIOR  ESTIMATION 

Figures  9  through  16  provide  the  pseudocode  used  for  our  posterior  estimation  algorithm. 


function  ( S ,  X,  6,  Y)  =  updateBackground(I,  S ,  X,  S,  Y,  M,  B ,  0,  /) 
for  all  regions  j,  pixels  p  £  Qj  do 
%  Stationary  component 


rk  =  Ei= i 


I  —  X  —  S.  *Y  —  D.  *  M 


-1  k,f,i 


(p)  k  =  1,2, ...  ,K 


{Skj(p)}k  -  A /"  (/x,  S);  (y,£)  =calcPost(S6(j),{rfc}fe) 


for  all  passes  i  do 
%  Speckle  noise 
rk  =  \l-  S-S.  *Y-D.  *M 


-I  k,f,i 


( P )  /c=  1,2,...,  if 


{^fc,/d(p)}fe  ~  AA  (/^,  S);  (y,£)  =calcPost(£x(j),{rfc}fc) 


%  Specular  noise 


rk  = 


I  -  S  -  X  -  D.*M 


kj,i 


(p)  = 


~A/'(y,£);  (y,£)  =calcPost(5Y  (j),{rfe}fc) 
end  for  %  (passes) 
end  for  %  (regions) 


(p),  r,;  =  [5ft{rM}fc  ;S{rfe,Jfc]i , 


%  Specular  noise  indicators 
rk,i  =  I— S  —  X  —  D.*M 

r  =  (Eti  nJ/JV,  [y]fc  =  (Eii’n,/,i(p))/JV 

q0  =  NormPDF(r;0,Ev/7V2),  gi  =  NormPDF(r;  [9ft{y}  ;  3{y}]T,  (£y  +  Er)/iV2) 


<*/(p)  ~  Bernoulli  (^-) 


end  function 


Figure  10.  Background  update 


function  (£>,  M)  =  updateForeground(/,  5,  X,  <5,  Y,  M,  B,  0,  /) 
for  all  passes  i  do 

%  Target  amplitudes 


rk  = 


I  —  S  —  X  —  5.  *Y 


{Mkj,i(p)}k 


>)  fc=l,2,...,Ff 

■A/>,E);  (m,S)  =calcPost(£M,{rfc}A 


%  Specular  noise  indicators 
r  =  ;S{rfc}fe]T,  [m]fc  =  Mk,f,i(p ) 

90  =  NormPDF(r;  0,  £v ) 

91  =  NormPDF(r;  [5ft  {y}  ;  3  {y}]T,  Ev  +  EM) 
Df,i(P)  ~  Bernoulli  (^-) 

end  for  %  (passes) 
end  function 


Figure  11.  Foreground  update 


function  0  =  updateHyperparameters(/,  S,  X ,  S ,  Y,  M,  B,  0) 

%  Region  classification 

(d,  q ,  z)  =  updateRegionClassifications(/,  S ,  X,  S,  Y,  M,  B,  0) 

%  BACKGROUND  covariance  updates 
for  all  regions  j  do 

%  Stationary  covariance 

sf  =  complex2real(S'fej/(p)),  V/;  s  =  concatenate({s/}j) 
(sS(j).a'S(i)>PS(i))  =  updateCovariance(CTS(j),  ssT ,  z(j)  *  F) 

%  Speckle  covariance 

Xfti  =  complex2real(XfcI/,»(p)),  V/,  i;  x  =  concatenate J 
(S-Y(j),  ax(j ),  px(j))  =  updateCovariance(erY  (j ),  px (j),  xxT ,  z(j)FN) 

%  Specular  covariance 

j7/,»  =  complex2real(Yfci/)»(p)),  Vf,i;  y  =  concatenate^/, j}^) 
{T,Y(j),a^  {j),py  (j))  =  updateCovariance(o-i  (j),p*  (j),yyT ,z(j)FN) 

end  for 


%  OTHER  covariance  matrices 
%  Target  covariance 

m/,i  =  complex2real(Mfcjji(p)),  V/,  i;  m  =  concatenate({m/y}^  J 
(SM,  ctm)  =  updateCovariance(erM ,  0,  mmT ,  NxNyFN) 

%  Additive  noise 

Vkj,i(p)  =  hj,i(p)  -  Skj{p)  -  XkJ:i(p)  -  Sf(p).  *  Yk>fti(p)~  Dfti(p).  *  MkJti(p) 
Vfti  =  complex2real(T4jii(p)),  V/,  i;  v  =  concatenate^?;/, 

(Y,v ,crv)  =  updateCovariance(cr'/ ,  0,  vvT ,  NxNyFN) 


%  INDICATOR  PROBABILITIES 

{7r/fiW}/.  =updateIndicatorProbabilities(D,e^QtiQi,e^mporai,7V)  %  Target 
| irj (p) |  =  updateIndicatorProbabilities((5,  expatial,  ^Temporal’  1)  %  Specular  Noise 


Q  =  l[d,q,{Xs(j),?:x(j)^(j)}j 


end  function 


Figure  12.  Hyper-parameters  update 


function  (d,  q,  z)  =  updateRegionClassifications(/,  S,  X,  S ,  Y,  M,  B ,  0) 
for  all  pixels  p  do 
for  all  classes  j  do 

Sf  =  complex2real(5'fci/(p)),  V/;  s  =  concatenate  ({s/}^) 

Xfti  =  complex2real(Xfc  V/,*;  x  =  concatenated:?/,^  J 

Ts  =  —trace  ([T,s (j)]~1ssT) ,  Ns  =  F 

Tx  =  -trace  ([Ex  (j)]~1xxT) ,  Nx  =  FN 

Wj(p)  =  exp  {Ts  -  Ns/2\og\T,s(j)\  -  KNS  log(27r)  +TX  -  Nx/2  log  |Sv(j)|  -  KNX  log(27r)  + 

end  for 
end  for 

if  useHMM  then 

Wj{p)  =  HMM({wj(p)}p  ,p) 

end  if 

%  Class  Assignment 
for  all  pixels  p  do 

d{p)  =  Multinomial (1, 

end  for 

%  Prior  probabilities 

for  all  classes  j  do 

z(j)  =  I  {Pi  dj(p)  =  1}  |  %  Number  of  pixels  with  class  j 

end  for 

q  =  Gamma(?+  c,  1  )/C 

end  function 


Figure  13.  Region  classification  update  pseudocode 


function  (S,  a.  p)  =  updateCovariance(<r,  p,  T ,  L) 
a2  =  Inv  —  Gamma  {aa  +  L,ba  +  trace(T)/2) 

£  =  scaleMat  ( cr2,  p) 

Pbest  =  Inv  —  WishartPDFCE ;  T  +  scaleMat(a2,  p)  *  ( df),df  +  N  +  2K  +  1) 

for  all  Metropolis-Hastings  repeats  do 
Pnew  —  Beta(ap,  &p) 

Pnew  =  Inv  —  WishartPDF{H ;  T  +  scaleMat(er2,  pnew)  *  ( df),df  +  N  +  2 K  +  1) 

if  Uniform (0, 1)  <  Pnew/Pbest  then 
P  —  Pnew  i  Pbest  —  Pnew 

end  if 
end  for 
end  function 

function  £  =  scaleMat(<r2,  p) 

r  =  2<t2  (pllT  +  (l-p)/*) 

s=  [  r  0  1 
^  [or 

end  function 

Figure  14.  Covariance  update  pseudocode 


function  n  =  updateIndicatorProbabilities(.B,  e spatial,  ^temporal,  L) 

for  l  =  1  to  L  do 

for  all  pixels  p,  frames  /  do 

^ high-)  (jPi  ^f,i)  ^  £ spatial  cind  W f  ^i(jp  ^  B  f >  £ temporal 

II 

S3 

^  c^iow  ->  else 

bhighi  ^ spatial  &nd  ~Wf,i  (jPi  l,i)  ^  £ temporal 

bfAp )  =  j 

blow  5  else. 

7T/,i(p)  =  Beta (a/,/(p)  +  bfj(p)  +  1  -  5/,i(p)) 

end  for 

end  for 

end  function 

Figure  15.  Indicator  probability  update  pseudocode 


function  (/x,  £)  =  calcPost(Sprior.,  (?’fc}^=1) 

A  _  y’  ( y  I 

—  ^jpri0r  y^-Jprior  \  ) 

r  =  complex2real({rfc}^=1) 

/x  =  At*;  5]  =  Sprjor ( I2k  A) 

end  function 

function  r  =  complex2real({rfe}^=1) 
r  =  [5RK}fc;3{rfc}fc]T 

end  function 

function  R  —  concatenate ({T*n}^=1) 

R=[ri  r2  ■■■  tat] 


end  function 


Figure  16.  Other  functions  used  in  this  pseudocode 


